home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The 640 MEG Shareware Studio 2
/
The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO
/
pascal
/
gaussj.zip
/
GAUSSJ.LIB
< prev
Wrap
Text File
|
1985-04-03
|
3KB
|
167 lines
{87}
procedure gaussj
(var b: ary2s; { square matrix of coefficients }
y: arys; { constant vector }
var coef: arys; { solution vector }
ncol: integer; { order of matrix }
var error: boolean); { true if matrix singular }
{ Gauss Jordan matrix inversion and solution }
{ B(n,n) coefficient matrix becomes inverse }
{ Y(n) original constant vector }
{ W(n,m) constant vector(s) become solution vector }
{ DETERM is the determinant }
{ ERROR=1 if singular }
{ INDEX(n,3) }
{ NV is number of constant vectors }
label 99;
var
w : array[1..maxc,1..maxc] of real;
index : array[1..maxc,1..3] of integer;
i,j,k,l,nv,
irow,icol,
n,l1 : integer;
determ,pivot,
hold,sum,t,
ab,big : real;
procedure swap(var a,b: real);
var hold : real;
begin { swap }
hold:=a;
a:=b;
b:=hold
end; { procedure swap }
procedure gausj2;
label 98;
var i,j,k,l,l1 : integer;
procedure gausj3;
var l : integer;
begin { procedure gausj3 }
{ interchange rows to put pivot on diagonal }
if irow<>icol then
begin
determ:=-determ;
for l:=1 to n do
swap(b[irow,l],b[icol,l]);
if nv>0 then
for l:=1 to nv do
swap(w[irow,l],w[icol,l])
end { if iroe<>icol }
end; { gausj3 }
begin { procedure gausj2 }
{ actual start of gaussj }
error:=false;
nv:=1; { single constant vector }
n:=ncol;
for i:=1 to n do
begin
w[i,1]:=y[i]; { copy constant vector }
index[i,3]:=0
end;
determ:=1.0;
for i:=1 to n do
begin
{ search for largest element }
big:=0.0;
for j:=1 to n do
begin
if index[j,3]<>1 then
begin
for k:=1 to n do
begin
if index[k,3]>1 then
begin
writeln('ERROR: matrix is singular');
error:=true;
goto 98 { abort }
end;
if index[k,3]<1 then
if abs(b[j,k])>big then
begin
irow:=j;
icol:=k;
big:=abs(b[j,k])
end
end { k-loop }
end
end; { j-loop }
index[icol,3]:=index[icol,3]+1;
index[i,1]:=irow;
index[i,2]:=icol;
gausj3; { further subdivision of gaussj }
{ divide pivot row by pivot column }
pivot:=b[icol,icol];
determ:=determ*pivot;
b[icol,icol]:=1.0;
for l:=1 to n do
b[icol,l]:=b[icol,l]/pivot;
if nv>0 then
for l:=1 to nv do
w[icol,l]:=w[icol,l]/pivot;
{ reduce nonpivot rows }
for l1:=1 to n do
begin
if l1<>icol then
begin
t:=b[l1,icol];
b[l1,icol]:=0.0;
for l:=1 to n do
b[l1,l]:=b[l1,l]-b[icol,l]*t;
if nv>0 then
for l:=1 to nv do
w[l1,l]:=w[l1,l]-w[icol,l]*t;
end { if l1<>icol }
end
end; { i-loop }
98:
end; { gausj2 }
begin { gaus-jordan main program }
gausj2; { first half of gaussj }
if error then goto 99;
{ interchange columns }
for i:=1 to n do
begin
l:=n-i+1;
if index[l,1]<>index[l,2] then
begin
irow:=index[l,1];
icol:=index[l,2];
for k:=1 to n do
swap(b[k,irow],b[k,icol])
end { if index }
end; { i-loop }
for k:=1 to n do
if index[k,3]<>1 then
begin
writeln(chr(7),'ERROR: matrix is singular');
error:=true;
goto 99 { abort }
end;
for i:=1 to n do
coef[i]:=w[i,1];
99:
end; { procedure gaussj }